# TODO: Add comment
# 
# Author: yaping
###############################################################################


for (e in commandArgs(TRUE)) {
	ta = strsplit(e,"=",fixed=TRUE)
	if(! is.na(ta[[1]][2])) {
		if(ta[[1]][1] == "wd"){
			wd<-ta[[1]][2]  ## directory to make plot file
			
		}
		if(ta[[1]][1] == "file"){
			file<-ta[[1]][2]
		}
		if(ta[[1]][1] == "anno_head"){
			anno_head<-ta[[1]][2]
		}
		if(ta[[1]][1] == "genome_size"){
			genome_size<-ta[[1]][2]
		}

	}
}

setwd(wd)
#file="NS157_MAR_cluster_merge500bp.fdr0.01.bed.annotate_detail.txt"
anno<-read.table(file,sep="\t",header=F)
#anno_head="NS157_MAR_cluster_merge500bp.fdr0.01.bed.annotate_detail.head.txt"
head_file<-read.table(anno_head,sep="\t",header=F)
head<-head_file[,1:(length(head_file[1,])-1)]
name<-head[,seq(1,length(head[1,]),by=3)]
value<-head[,seq(2,length(head[1,]),by=3)]
feature_total<-head[,seq(3,length(head[1,]),by=3)]

feature_overlap_count<-colSums(anno[,seq(5,length(anno[1,]),by=2)],na.rm=T)
feature_overlap_count<-feature_overlap_count[1:(length(feature_overlap_count)-1)]
feature_overlap_percentage<-feature_overlap_count/feature_total
count<-c()
for(i in seq(5,length(anno[1,]),by=2)){
	count<-c(count,dim(anno[anno[,i]!=0,])[1])
}
count<-count[1:(length(count)-1)]
percentage<-count/length(anno[,1])
overlap_base<-colSums((anno[,3]-anno[,2])*anno[,seq(6,length(anno[1,]),by=2)],na.rm=T)
ndr_length<-sum(anno[,3]-anno[,2], na.rm=T)
#p<-phyper(overlap_base, t(value), genome_size-t(value), ndr_length, log.p=T, lower.tail=F)
#p<--p/log(10)

#sum<-cbind(t(name),count,percentage,feature_overlap_count,t(feature_overlap_percentage),p)
sum<-cbind(t(name),count,percentage,feature_overlap_count,t(feature_overlap_percentage))

write.table(sum,paste(file,".annotation.summary.txt",sep=""),quote=F, sep="\t",row.names=F,col.names=F)
